clear all;close all;clc
fs = 44.1e3;
f1 = 200;
f2 = 1000;
f3 = 12000;

figure('position',[80 80 700 300]);

xv =[5 9 17 21] ; % horizontal 1-25 locations
yv = 1:5:25; % vertical 1-50 locations
x_loc = [-80:15:-65 -55 -45:5:45 55 65:15:80]; %Azimuths
y_loc = -45:5.625:230.625; %Elevations

for One1_Two2 = 2:-1:1 % 2 plot monaural HRTFs   1; % plot binaural HRTFs
    
    folder = 'standard_hrir_database/';
    SN = {'003';'008';'009';'010';'011';'012';'015';'017';'018';'019';'020';'021';'027';...
        '028';'033';'040';'044';'048';'051';'058';'059';'060';'061';'065';'119';'124';...
        '126';'127';'131';'133';'134';'135';'137';'147';'148';'152';'153';'154';'155';'156'};
    X_ILD = [];X_ITD = [];
    d = [];
    c_row = 0;
    
    NC = length(xv)*length(yv);
    c = 0;
    for azn=xv
        c_row=c_row+1;
        c_col = 0;
        for eln=yv
            c_col=c_col+1;
            c = c+1;
            if One1_Two2 == 1
                Bianural_vectors{c} = [];
                Bianural_vectors{c} = [];
            else
                Monaural_vectors{c_row,c_col} = [];
                Monaural_vectors{c_row,c_col} = [];
            end
        end
    end
    
    for isub = 1:1
        isub
        c = 0;
        c_row = 0;
        eval(['load ' folder 'subject_' SN{isub} '/hrir_final'])
        for azn=xv
            azn
            c_col = 0;
            c_row=c_row+1;
            temp_vect = [];
            for eln=yv
                c_col=c_col+1;
                h_l = squeeze(hrir_l(azn,eln,:));
                h_r = squeeze(hrir_r(azn,eln,:));
                Nfft = length(h_l);
                fsig = (1:Nfft)./Nfft.*fs;
                N1 = find(fsig>f1);N1 = N1(1);
                N2 = find(fsig<=f2);N2 = N2(end);
                N3 = find(fsig<=f3);N3 = N3(end);
                yl = abs(fft(h_l,Nfft))';
                yr = abs(fft(h_r,Nfft))';
                phasel = angle(fft(h_l,Nfft))';
                phaser = angle(fft(h_r,Nfft))';
                if One1_Two2==1
                    phrtf =  [(phaser(N1:N2)-phasel(N1:N2)) yr(N2+1:N3)-yl(N2+1:N3)];
                else
                    phrtf =  yl(N1:N3);
                end
                temp_vect=[temp_vect
                    phrtf];
            end % elevation finished
            subplot(2,5,c_row+ (3-One1_Two2-1)*4)
            if One1_Two2==1 % binaural
                temp_vect=f_HRTF_wrap(temp_vect);
                m = 0;
                for kk = c+1:c+length(yv)
                    m = m+1;
                    Bianural_vectors{kk}=[Bianural_vectors{kk}
                        temp_vect(m,:)];
                end
                c = c+length(yv);
                
                for elevation = 1:length(yv)
                    semilogx(fsig(N1:N3),temp_vect(elevation,:)); hold on;
                end
                plot([1000 1000],[-5 5],'k:');
                title(['Hori = ' num2str(x_loc(azn)) ' deg'])
                axis([150 15000 -5 5])
                set(gca,'xtick',[200 1000 10000],'xticklabel',[200 1000 10000],'ytick',-5:2.5:5)
                box off
            else % monaural
                for elevation = 1:length(yv)
                    semilogx(fsig(N1:N3),20*log10(temp_vect(elevation,:)),'.-'); hold on;
                end
                title(['Hori = ' num2str(x_loc(azn)) ' deg'])
                axis([150 15000 -30 15])
                set(gca,'xtick',[200 1000 10000],'xticklabel',[200 1000 10000])
                box off
                
            end
        end
        pause(0.1)
    end % subject
end
for ileg = 1:length(yv)
    yleg{ileg} = ['Vert: ' num2str(round(y_loc(yv(ileg)))) ' deg'];
end

legend(yleg)







